{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"provenance":[{"file_id":"1PC10pGNfjXTJz0p3RrCivWC2bwxcprOj","timestamp":1715631212751}],"collapsed_sections":["p0jrcMSbKdix","CVwOxESDKpS_"],"authorship_tag":"ABX9TyMArFvtbiV9yqbn9YMIUi7o"},"kernelspec":{"name":"python3","display_name":"Python 3"},"language_info":{"name":"python"}},"cells":[{"cell_type":"code","source":["# Install packages\n","import pandas as pd\n","import numpy as np"],"metadata":{"id":"YqgTZcJn3bP2"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Mount my Google Drive\n","from google.colab import drive\n","drive.mount('/content/drive')"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"kXDXoLGtKEk9","executionInfo":{"status":"ok","timestamp":1715807831175,"user_tz":420,"elapsed":1204,"user":{"displayName":"Amanda Smith","userId":"00495916758981390713"}},"outputId":"704a7748-170b-47e2-d495-ae16197093fa"},"execution_count":null,"outputs":[{"output_type":"stream","name":"stdout","text":["Drive already mounted at /content/drive; to attempt to forcibly remount, call drive.mount(\"/content/drive\", force_remount=True).\n"]}]},{"cell_type":"markdown","source":["# Functions & Classes"],"metadata":{"id":"SBxxCw7HO648"}},{"cell_type":"code","source":["# Round off (for starting sizing)\n","def round_to_nearest(number, round_to):\n","    \"\"\"Rounds a number to the nearest multiple. Default 10\"\"\"\n","    return round(number / round_to) * round_to"],"metadata":{"id":"T-sZ3t-6O-hG"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Define battery class\n","class Battery:\n","    def __init__(self, capacity, efficiency, charge_rate, discharge_rate, init_charge=0):\n","        \"\"\"Initialize battery parameters.\"\"\"\n","        self.capacity = capacity * 0.95   # Total energy storage capacity (kWh)\n","        # with depth of discharge limited to 80% for battery longevity\n","        self.efficiency = efficiency  # Round-trip efficiency (0.0 - 1.0)\n","        self.charge = init_charge     # Initial state of charge (kWh)\n","        self.max_charge_rate = charge_rate * capacity  # Maximum charge rate (kW)\n","        self.max_discharge_rate = discharge_rate * capacity  # Maximum discharge rate (kW)\n","\n","    def update_charge(self, demand, pv_output):\n","        \"\"\"\n","        Update the battery charge based on demand and PV output.\n","\n","        \"\"\"\n","\n","        # Calculate net energy balance\n","        net_energy = pv_output - demand\n","\n","        if net_energy > 0:  # Surplus energy\n","            max_charge = min(net_energy, self.max_charge_rate, self.capacity - self.charge)\n","            self.charge += max_charge * self.efficiency\n","        else:  # Energy deficit\n","            max_discharge = min(-net_energy, self.max_discharge_rate, self.charge)\n","            self.charge -= max_discharge / self.efficiency\n","\n","        # Ensure charge stays within bounds\n","        self.charge = max(0, min(self.charge, self.capacity))"],"metadata":{"id":"zfzUS3ceQkHv"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["# Population sizes"],"metadata":{"id":"wTuVG3NBl_uT"}},{"cell_type":"code","source":["# VARIABLES\n","\n","population_size = [25, 50, 500, 1000, 1500, 3000]\n","num_person = 8"],"metadata":{"id":"QqDLfn-ymDTi"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["# Electric Load"],"metadata":{"id":"p0jrcMSbKdix"}},{"cell_type":"code","source":["# ASSUMPTIONS\n","\n","elec_use_daily_person = 1.5 # (see lit notes)\n","household_size = 8          # 8.3 average for Senegal (see lit notes)"],"metadata":{"id":"PwzUTDOAH9_2"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Path to load profile in Google Drive\n","load_path = '/content/drive/My Drive/1 - projects/Senegal USAID/data/Load/Huld 2017 Fig. 4 Low nighttime consumption profile.csv'"],"metadata":{"id":"6JjW-xZL15Ez"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Import load data from csv of extracted profile\n","data = np.genfromtxt(load_path, delimiter=',')\n","\n","x = data[:, 0]  # First column (index 0) is x (hr)\n","y = data[:, 1]  # Second column (index 1) is y (kW)"],"metadata":{"id":"RlqYv6OgN7JA"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Convert load data to a dataframe\n","load_raw = pd.DataFrame(data, columns=['x', 'y'])\n","\n","# Generate integer x values from 1 to 24 (for hr of the day)\n","# and convert to float for consistency with df column of x values\n","all_x_values = pd.Series(range(1, 25), name='x').astype(float)\n","\n","# Merge the given x values with all_x_values\n","merged_df = pd.merge(load_raw, all_x_values.to_frame(), on='x', how='outer')\n","\n","# Sort to get x-values in sequential order\n","merged_df = merged_df.sort_values('x')\n","\n","# Perform linear interpolation to fill NaN values for y\n","merged_df['y'] = merged_df['y'].interpolate(method='linear')\n","\n","# Extract the power consumption number, y, for integer x values from 1 to 24.\n","merged_df = merged_df[merged_df['x'].isin(all_x_values)]\n","\n","# Move 24 (midnight) value to beginning of time series\n","# and re-index\n","# and set the x-value to 0\n","midnight_row = merged_df.iloc[-1:]\n","other_rows = merged_df.iloc[:-1]\n","loadprofile_df = pd.concat([midnight_row, other_rows])\n","loadprofile_df = loadprofile_df.reset_index(drop=True)\n","loadprofile_df.at[0,'x'] = 0.0\n","\n","# Convert column x to integers\n","loadprofile_df['x'] = loadprofile_df['x'].astype(int)"],"metadata":{"id":"iTDpqhWPymik"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Verify that it is scaled to 300 kWh/day consumption per Huld (2017).\n","original_scale = sum(loadprofile_df['y'])\n","\n","original_scale"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"Hgg9dUmPzaMq","executionInfo":{"status":"ok","timestamp":1715807831175,"user_tz":420,"elapsed":3,"user":{"displayName":"Amanda Smith","userId":"00495916758981390713"}},"outputId":"1821df0f-8370-4ec7-d7f3-197e7073da27"},"execution_count":null,"outputs":[{"output_type":"execute_result","data":{"text/plain":["299.7066244563945"]},"metadata":{},"execution_count":10}]},{"cell_type":"code","source":["# Rescale to different population sizes\n","scaling_factor = elec_use_daily_person * num_person / original_scale\n","loadprofile_df['y'] *= scaling_factor"],"metadata":{"id":"X7wTHh0FIMvb"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Make into a yearly load profile by repeating 365x.\n","yearly_loadprofile = np.tile(loadprofile_df.values, (365, 1))\n","yearly_loadprofile_df = pd.DataFrame(yearly_loadprofile)\n","\n","# Name columns\n","# assuming power value is constant through the hour, so that demand = kW * 1h\n","yearly_loadprofile_df.columns = ['Hr of day','Demand (kWh)']"],"metadata":{"id":"ZXSgClae__Mt"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["# PV System"],"metadata":{"id":"CVwOxESDKpS_"}},{"cell_type":"code","source":["# ASSUMPTIONS\n","\n","panel_eff = 0.19    # PVWatts default for Crystalline Silicon\n","tilt = 15           # close to latitude\n","orientation = 180   # S"],"metadata":{"id":"0tgSOe41K166"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Size of PV array simulated in PVWatts\n","sim_size_kWp = 10\n","\n","# Path to PV output profile in Google Drive\n","PVdata_path = '/content/drive/My Drive/1 - projects/Senegal USAID/data/PV/pvwatts_hourly_Dakar_10kW_groundmount.csv'"],"metadata":{"id":"IMdUnOS7BIh6"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Import load data from csv of extracted profile\n","PVoutput_df = pd.read_csv(PVdata_path, skiprows=31)"],"metadata":{"id":"EmwHkjmB3wYG"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Initial guess for PV sizing\n","daily_demand = sum(yearly_loadprofile_df['Demand (kWh)'])/365\n","pv_sizing_start = daily_demand/4  # kWp\n","\n","# Allow size increments of 1 kW below 25 kW, 5 kW between 25-100 kW,\n","# and 10 kW above 100 kW.\n","if pv_sizing_start < 25:\n","    size_kWp = round_to_nearest(pv_sizing_start,1)\n","elif pv_sizing_start < 100:\n","    size_kWp = round_to_nearest(pv_sizing_start,5)\n","else:\n","    size_kWp = round_to_nearest(pv_sizing_start,10)\n","\n","# NOTE initial guess must be lower than the appropriate PV array size, or\n","# the sizing algorithm should be modified to avoid unnecessary oversizing.\n","\n","size_kWp"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"zWQKHWJ6kgpR","executionInfo":{"status":"ok","timestamp":1715807831520,"user_tz":420,"elapsed":6,"user":{"displayName":"Amanda Smith","userId":"00495916758981390713"}},"outputId":"65c19452-1f23-4d0a-e306-9273630e09d0"},"execution_count":null,"outputs":[{"output_type":"execute_result","data":{"text/plain":["3"]},"metadata":{},"execution_count":16}]},{"cell_type":"code","source":["# Size PV array\n","\n","for i in range(200):\n","\n","    # Scale output according to PV system size\n","    PVoutput_df['AC System Output (W)'] = PVoutput_df['AC System Output (W)'] * size_kWp / sim_size_kWp\n","    PVoutput_df['DC Array Output (W)'] = PVoutput_df['DC Array Output (W)'] * size_kWp / sim_size_kWp\n","\n","    # Convert AC system output (W) values to kWh values for the previous hour\n","    # by converting W to kW and assuming power provided is constant\n","    # through the previous hour, so that PV output = kW * 1h\n","    PVoutput_df['PV Output (kWh)'] = PVoutput_df['AC System Output (W)'] / 1000\n","\n","    # Check whether the PV output exceeds annual demand by 15%\n","    # (This value was found through trial & error and will be consistent\n","    # for the given load profile.)\n","    sufficient_pv = sum(PVoutput_df['PV Output (kWh)']) > (sum(yearly_loadprofile_df['Demand (kWh)']) * 1.15)\n","\n","    if sufficient_pv:\n","        print(f\"Sufficient PV size at iteration {i} is: {round(size_kWp,2)} kWh\")\n","        break\n","\n","    else:\n","        # Allow size increments of 1 kW below 25 kW, 5 kW between 25-100 kW,\n","        # and 10 kW above 100 kW.\n","        if size_kWp < 25:\n","            size_kWp += 1\n","        elif size_kWp < 100:\n","            size_kWp += 5\n","        else:\n","            size_kWp += 10\n","\n","    print(i)\n","    print(size_kWp)\n","\n","else:\n","    print(f\"Did not find sufficient PV size in {i} iterations.\")"],"metadata":{"id":"Dj-C9yPZW74H","executionInfo":{"status":"ok","timestamp":1715807831520,"user_tz":420,"elapsed":5,"user":{"displayName":"Amanda Smith","userId":"00495916758981390713"}},"colab":{"base_uri":"https://localhost:8080/"},"outputId":"680b5337-e3c8-4098-968e-3191936974bd"},"execution_count":null,"outputs":[{"output_type":"stream","name":"stdout","text":["1.0946152246575638\n","0\n","4\n","0.4378460898630259\n","1\n","5\n","0.21892304493151296\n","2\n","6\n","0.1313538269589074\n","3\n","7\n","0.09194767887123508\n","4\n","8\n","0.07355814309698855\n","5\n","9\n","0.06620232878728957\n","6\n","10\n","0.06620232878728957\n","7\n","11\n","0.07282256166601848\n","8\n","12\n","0.08738707399922213\n","9\n","13\n","0.11360319619898877\n","10\n","14\n","0.15904447467858424\n","11\n","15\n","0.23856671201787635\n","12\n","16\n","0.38170673922860004\n","13\n","17\n","0.6489014566886252\n","14\n","18\n","1.168022622039525\n","Sufficient PV size at iteration 15 is: 18 kWh\n"]}]},{"cell_type":"code","source":["# Calculate array area\n","PV_area = size_kWp / panel_eff # kW * (1 m2/kW under test conditions) * kW/kW\n","PV_area # (m2)"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"fDUKA-17LgtW","executionInfo":{"status":"ok","timestamp":1715807831520,"user_tz":420,"elapsed":5,"user":{"displayName":"Amanda Smith","userId":"00495916758981390713"}},"outputId":"3faab83c-6366-4ea4-d939-112dc4ab62f9"},"execution_count":null,"outputs":[{"output_type":"execute_result","data":{"text/plain":["94.73684210526315"]},"metadata":{},"execution_count":18}]},{"cell_type":"markdown","source":["# Combine working data into 1 dataframe and compare PV output with demand"],"metadata":{"id":"10wChh4fDAeN"}},{"cell_type":"code","source":["# Filter load dataframe for demand column\n","yearly_loadprofile_df_filtered = yearly_loadprofile_df[['Demand (kWh)']]\n","\n","# Filter PV output dataframe for time-related and PV output columns\n","PVoutput_df_filtered = PVoutput_df[['Month','Day','Hour','PV Output (kWh)']]\n","\n","# Concatenate into 1 df\n","df = pd.concat([PVoutput_df_filtered, yearly_loadprofile_df_filtered], axis=1)"],"metadata":{"id":"Kh1q56W_DG8k"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Compare total PV output to total demand.\n","PVoutput_to_demand_annual = sum(df['PV Output (kWh)'])/sum(df['Demand (kWh)'])\n","PVoutput_to_demand_annual"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"dCOBwUBALyIs","executionInfo":{"status":"ok","timestamp":1715807831520,"user_tz":420,"elapsed":4,"user":{"displayName":"Amanda Smith","userId":"00495916758981390713"}},"outputId":"714d70da-1d4b-4694-927f-4331239e15c1"},"execution_count":null,"outputs":[{"output_type":"stream","name":"stdout","text":["5115.939084532981\n","4379.999999999882\n"]},{"output_type":"execute_result","data":{"text/plain":["1.168022622039525"]},"metadata":{},"execution_count":20}]},{"cell_type":"markdown","source":["# Battery and check reliability"],"metadata":{"id":"aPs432oDKpyM"}},{"cell_type":"code","source":["# ASSUMPTIONS\n","\n","efficiency = 0.95     # round-trip efficiency\n","\n","charge_rate = 0.8     # (capacity/hr)\n","discharge_rate = 1.0  # (capacity/hr)\n","\n","reliability_required = 0.95    # portion of hours that demand should be met by PV + batt"],"metadata":{"id":"Ymqcg8YfTT8S"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Initial guess for battery sizing\n","battery_sizing_start = daily_demand/3\n","\n","# Allow size increments of 5 kWh below 100 kWh, 10 kWh above.\n","if battery_sizing_start < 100:\n","  capacity = round_to_nearest(battery_sizing_start,5)         # (kWh)\n","else:\n","  capacity = round_to_nearest(battery_sizing_start,10)         # (kWh)\n","\n","# NOTE initial guess must be lower than the appropriate battery size, or\n","# the sizing algorithm should be modified to avoid unnecessary oversizing."],"metadata":{"id":"1lGA_mNkaCLN"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Size battery (warning--slow!)\n","for i in range(200):\n","\n","  # Instantiate battery\n","  battery = Battery(capacity=capacity, efficiency=efficiency, charge_rate=charge_rate, discharge_rate=discharge_rate)\n","\n","  # Simulate battery operation over each row (hourly timestep)\n","  for index, row in df.iterrows():\n","      battery.update_charge(row['Demand (kWh)'], row['PV Output (kWh)'])\n","      df.at[index, 'Charge (kWh)'] = battery.charge\n","\n","  # Calculate battery output\n","  df['Previous Charge (kWh)'] = df['Charge (kWh)'].shift(1).fillna(0)\n","  df['Battery Output (kWh)'] = df['Previous Charge (kWh)'] - df['Charge (kWh)']\n","\n","  # Identify hours where PV + battery combined do not meet the demand\n","  df['Insufficiency'] = df['Demand (kWh)'] > (df['PV Output (kWh)'] + df['Battery Output (kWh)'])\n","\n","  # Check whether it meets reliability criteria\n","  # with tolerance so that reliability rounds off to 95%\n","  sufficient = sum(df['Insufficiency']) <= (8760 * (1-(reliability_required-0.005)))\n","\n","  # Quantify reliability\n","  hours_met = (1-sum(df['Insufficiency']))/8760\n","  sim_reliability = (1 - sum(df['Insufficiency'])/8760)*100\n","\n","  if i==0:\n","      print(f\"Initial sizing run with battery size {round(capacity,0)} yields reliability of {round(sim_reliability,2)}.\")\n","\n","  if sufficient:\n","    print(f\"Sufficient battery size at iteration {i} is: {round(capacity,2)} kWh\")\n","    break\n","\n","  else:\n","    # Allow size increments of 5 kWh below 100 kWh, 10 kWh above.\n","    if capacity < 100:\n","      capacity += 5\n","    else:\n","      capacity += 10\n","\n","else:\n","  print(f\"Did not find sufficient battery size in {i} iterations.\")\n","  print(f\"Current battery size is {round(capacity,0)} kWh and reliability is {round(sim_reliability,2)}.\")"],"metadata":{"id":"F45pRnOjPWiH","executionInfo":{"status":"ok","timestamp":1715807837414,"user_tz":420,"elapsed":5897,"user":{"displayName":"Amanda Smith","userId":"00495916758981390713"}},"colab":{"base_uri":"https://localhost:8080/"},"outputId":"e1beb85a-b206-460e-fc28-bbec976b7442"},"execution_count":null,"outputs":[{"output_type":"stream","name":"stdout","text":["Initial sizing run with battery size 5 yields reliability of 80.24.\n","Sufficient battery size at iteration 1 is: 10 kWh\n","95.4337899543379\n"]}]},{"cell_type":"markdown","source":["# Cost Estimate"],"metadata":{"id":"SpqgrioGKqFq"}},{"cell_type":"code","source":["# ASSUMPTIONS\n","\n","PV_cost_per_kWp = (950 + 1150)        # USD, from Szabo et al. 2021, including modules + installation and BOS\n","battery_cost_per_kWh = 400            # USD, from Szabo et al. 2021, for Li-ion\n","\n","# NOTE\n","# O&M costs (for battery and panel maintenance) are ignored here!"],"metadata":{"id":"NYBwZpfJ3N4r"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Estimate PV system cost\n","PV_cost = PV_cost_per_kWp * size_kWp\n","PV_cost"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"9yR0xWBD5zOC","executionInfo":{"status":"ok","timestamp":1715807837414,"user_tz":420,"elapsed":21,"user":{"displayName":"Amanda Smith","userId":"00495916758981390713"}},"outputId":"7a8a3f06-dba4-478e-fdec-b2f541fc92dd"},"execution_count":null,"outputs":[{"output_type":"execute_result","data":{"text/plain":["37800"]},"metadata":{},"execution_count":25}]},{"cell_type":"code","source":["# Estimate battery cost\n","battery_cost = battery_cost_per_kWh * capacity\n","battery_cost"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"HKq1Bby06H51","executionInfo":{"status":"ok","timestamp":1715807837414,"user_tz":420,"elapsed":19,"user":{"displayName":"Amanda Smith","userId":"00495916758981390713"}},"outputId":"f8dbaec2-aa5e-4762-defb-37c371d753b8"},"execution_count":null,"outputs":[{"output_type":"execute_result","data":{"text/plain":["4000"]},"metadata":{},"execution_count":26}]},{"cell_type":"code","source":["# Estimate total system cost\n","total_cost = PV_cost + battery_cost\n","total_cost"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"0IG0uJ0b6ICf","executionInfo":{"status":"ok","timestamp":1715807837694,"user_tz":420,"elapsed":298,"user":{"displayName":"Amanda Smith","userId":"00495916758981390713"}},"outputId":"496cd835-6418-448b-d59a-11d793e19d2e"},"execution_count":null,"outputs":[{"output_type":"execute_result","data":{"text/plain":["41800"]},"metadata":{},"execution_count":27}]},{"cell_type":"markdown","source":["# Output file"],"metadata":{"id":"O-Of6An97r1S"}},{"cell_type":"code","source":["# Create output df\n","output_data = {\n","    'People': [num_person],\n","    'PV capacity (kWp)': [size_kWp],\n","    'PV panel area (m2)': [round(PV_area, 2)],\n","    'PV system cost ($)': [PV_cost],\n","    'Battery capacity (kWh)': [capacity],\n","    'Battery cost ($)': [battery_cost],\n","    'Total cost ($)': [total_cost]\n","}\n","output_df = pd.DataFrame(data=output_data)"],"metadata":{"id":"1wYfFM0o74LQ"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Create output file\n","output_file_name = f\"results_ground_tilt{tilt}deg_orientation{orientation}_{num_person}people.csv\"\n","output_file_path = \"/content/drive/My Drive/1 - projects/Senegal USAID/\"\n","output_df.to_csv(output_file_path + output_file_name, index=False)  # exclude row numbers"],"metadata":{"id":"nNQDw7zoSr2t"},"execution_count":null,"outputs":[]}]}